% flagImgExport = 0;

if ~exist('FTFITManual','var')
    FTFITManual = 1;
end
if FTFITManual
    paraSet = {'../', 12, '8nlps'};
    paraSetCellNum = 370;
end

fprintf('measurement = '); disp(paraSet{1});
fprintf('numbers = '); fprintf('%d  ',paraSet{2}); fprintf('\n');
fprintf('flowRate = '); disp(paraSet{3});
dataExtension = 'data_export';

measNum = paraSet{2};
flowRate = paraSet{3};
if ~exist('noOfFCs','var')
    noOfFCs = 10;
end

if FTFITManual
    flagImgShow = 1;
    flagImgExport = 1;
    flagInterrupt = 0;
else
    % flagImgExport = 1;
    flagImgShow = 0;
    flagInterrupt = 0;
end

filedir = [paraSet{1} filesep 'Online' filesep];
idx = find(filedir == filesep);
measName = filedir(idx(end-3)+1:idx(end-2)-1);
sampleName = filedir(idx(end-2)+1:idx(end-1)-1);
pathMeas = filedir(1:idx(end-2)-1);
if FTFITManual
    rootdir = [pathMeas filesep 'FTFitData' filesep 'CellsExport' filesep];
else
    rootdir = [pathMeas filesep 'FTFitData' filesep sprintf(['%s' '_' '%s' '_M%i-%i'], sampleName, flowRate, measNum(1), measNum(end)) filesep];
end

%% load data
load([filedir sprintf('M%i_%s', measNum(1)), dataExtension]);
sampleData = sample;
if FTFITManual
    mask = paraSetCellNum == [sampleData.measData.cellData.number];
    sampleData.measData.cellData(~mask) = [];
end
for n = 2:length(measNum)
    load([filedir sprintf('M%i_%s', measNum(n)), dataExtension]);
    sampleData(n) = sample;
end
clear sample;

if ~exist(rootdir,'dir')
    mkdir(rootdir);
end
if FTFITManual
    rootdir = sprintf([rootdir '%s' '_' '%s' '_'], sampleName, flowRate);
end

globalCount = 0;
for idxMeas = 1:length(measNum)
    globalCount = globalCount + length(sampleData(idxMeas).measData.cellData);
end

%% FT & fitting
subplott = @(m,n,p) subtightplot (m, n, p, [0.008 0.06], [0.09 0.10], [0.12 0.025]);

% output
defoTau1Inlet = zeros(1,globalCount);
defoTau1Outlet = zeros(1,globalCount);

defoTau2Inlet = zeros(1,globalCount);
defoTau2Outlet = zeros(1,globalCount);

taus1Circularity = zeros(1,globalCount);
taus2Circularity = zeros(1,globalCount);
rSqrsCircularityTau1 = zeros(1,globalCount);
rSqrsCircularityTau2 = zeros(1,globalCount);

tsMinDefo = zeros(1,globalCount);
tsPass = zeros(1,globalCount);

measNums = zeros(1,globalCount);
cellNums = zeros(1,globalCount);

divCellNumPtr = globalCount/20;
fprintf([repmat('-',1,20) '\n']);

fitFig = 100;
if ~ishandle(fitFig)
    figure(fitFig);
end
figPos = get(fitFig, 'Position'); set(fitFig, 'Position', [figPos(1), figPos(2), 570, 570]);
for k = 7:8
    if ~ishandle(k)
        figure(k); clf;
    end
end
phi = linspace(0,2*pi,101); phi(end) = []; phi = phi(:);
globalCount = 0;
inletGateHalfWidth = 3; % um
% measurement
channelPosMat = NaN(length(measNum),2);
for idxMeas = 1:length(measNum)
    channelPos = sampleData(idxMeas).channelPositions.um;
    channelPosMat(idxMeas,:) = channelPos;
    % cell
    for idxCell = 1:length(sampleData(idxMeas).measData.cellData)
        globalCount = globalCount + 1;
        fprintf(repmat('-',1,round(globalCount/divCellNumPtr) - round((globalCount-1)/divCellNumPtr)));
        cntrIdx = sampleData(idxMeas).measData.cellData(idxCell).data.contourIndex;
        circFCEven = zeros(1,length(cntrIdx));
        circFCOdd = zeros(1,length(cntrIdx));
        inertRatFCEven = zeros(1,length(cntrIdx));
        inertRatFCOdd = zeros(1,length(cntrIdx));
        % frame
        for idxFrame = 1:length(sampleData(idxMeas).measData.cellData(idxCell).data.time)
            cntr = sampleData(idxMeas).measData.cellData(idxCell).data.contourData(cntrIdx(idxFrame,1):cntrIdx(idxFrame,1)+cntrIdx(idxFrame,2)-1,:);
            if length(cntr) <= 2
                continue
            end
            % FT
            [FCa, FCp, ~, ak, bk] = funcCellContourFCManual(cntr(:,1), cntr(:,2), noOfFCs);
            % reconstruction
            FC = 1:2:noOfFCs-1;
            rFCOdd = repmat(ak(0+1),length(phi),1)...
                + sum( repmat(ak(FC+1),length(phi),1) .* cos(repmat(FC,length(phi),1).*repmat(phi,1,length(FC))) ,2)...
                + sum( repmat(bk(FC+1),length(phi),1) .* sin(repmat(FC,length(phi),1).*repmat(phi,1,length(FC))) ,2);
            
            FC = 2:2:noOfFCs-1;
            rFCEven = repmat(ak(0+1),length(phi),1)...
                + sum( repmat(ak(FC+1),length(phi),1) .* cos(repmat(FC,length(phi),1).*repmat(phi,1,length(FC))) ,2)...
                + sum( repmat(bk(FC+1),length(phi),1) .* sin(repmat(FC,length(phi),1).*repmat(phi,1,length(FC))) ,2);
            
            % calc circularity & inertia ratio
            circFCEven(idxFrame) = calcCircularity(phi, rFCEven);
            circFCOdd(idxFrame) = calcCircularity(phi, rFCOdd);
            inertRatFCEven(idxFrame) = calcInertiaRatio(phi, rFCEven);
            inertRatFCOdd(idxFrame) = calcInertiaRatio(phi, rFCOdd);
            
        end
        % mask for inside channel
        xPos = sampleData(idxMeas).measData.cellData(idxCell).data.xPosAcc;
        maskInsideChannel = (channelPos(1) <= xPos) & (xPos <= channelPos(2));
        [uniqueXPosAcc, idxA, ~] = unique(sampleData(idxMeas).measData.cellData(idxCell).data.xPosAcc);
        timeInsideChannel = interp1(uniqueXPosAcc, sampleData(idxMeas).measData.cellData(idxCell).data.time(idxA), channelPos);
        [~, idxA, ~] = unique(sampleData(idxMeas).measData.cellData(idxCell).data.time);
        maskInsideChannel = maskInsideChannel & any(idxA == (1:length(xPos)),1);
        time = sampleData(idxMeas).measData.cellData(idxCell).data.time - timeInsideChannel(1);
        tFine = linspace(0,timeInsideChannel(end)-timeInsideChannel(1),100);
        
        defoFCEven = 1-circFCEven(maskInsideChannel);
        % calc max deformation at channels inlet
        defoTau1Inlet(globalCount) = max(defoFCEven(1:3));
        try
            [para1, yFitCircEven, rSqr1] = expFit(time(maskInsideChannel),defoFCEven,-1,1,0,tFine);
        catch ME
            para1 = [mean(1-circFCEven(maskInsideChannel)) 0 -1];
            yFitCircEven = para1(2)*exp(-tFine/para1(3)) + para1(1);
            rSqr1 = 0;
        end
        % calc deformation at channels outlet = steady state of fit
        defoTau1Outlet(globalCount) = para1(1);
        
        try
            [para2, yFitCircOdd, rSqr2, defo2Inlet] = expFit_tau2(time(maskInsideChannel),1-circFCOdd(maskInsideChannel),+1,tFine);
        catch ME
            para2 = [mean(1-circFCOdd(maskInsideChannel)) 0 -1];
            yFitCircOdd = para2(2)*exp(-tFine/para2(3)) + para2(1);
            rSqr2 = 0;
            defo2Inlet = 0;
        end
        defoTau2Inlet(globalCount) = defo2Inlet;
        defoTau2Outlet(globalCount) = para2(1); % steady state = exp fit value for t -> infinity
        yFitCircAll = yFitCircEven + yFitCircOdd;
        
        tau1 = para1(3); tau2 = para2(3);
        taus1Circularity(globalCount) = tau1;
        taus2Circularity(globalCount) = tau2;
        rSqrsCircularityTau1(globalCount) = rSqr1;
        rSqrsCircularityTau2(globalCount) = rSqr2;
        
        tsPass(globalCount) = timeInsideChannel(end) - timeInsideChannel(1);
        tsMinDefo(globalCount) = log( -para1(2)/para2(2)*para2(3)/para1(3) ) * (para1(3)*para2(3))/(para2(3)-para1(3));
        measNums(globalCount) = measNum(idxMeas);
        cellNums(globalCount) = sampleData(idxMeas).measData.cellData(idxCell).number;
        
        if flagImgShow || flagImgExport
            set(0,'CurrentFigure',fitFig); clf;
            subplott(3,1,1);
                plot(time(maskInsideChannel),1-sampleData(idxMeas).measData.cellData(idxCell).data.circularity(maskInsideChannel),'.-', time(maskInsideChannel),1-circFCEven(maskInsideChannel)+1-circFCOdd(maskInsideChannel),'.-', tFine,yFitCircAll);
                title({sprintf('MeasNum = %i, cell.number = %i: start = %i, end = %i, length = %i', measNum(idxMeas), sampleData(idxMeas).measData.cellData(idxCell).number, sampleData(idxMeas).measData.cellData(idxCell).videoPointers(1), sampleData(idxMeas).measData.cellData(idxCell).videoPointers(end), sampleData(idxMeas).measData.cellData(idxCell).videoPointers(end)-sampleData(idxMeas).measData.cellData(idxCell).videoPointers(1)+1),...
                    ['\tau_1 = ' sprintf('%g ms, ',tau1) '\tau_2 = ' sprintf('%g ms, ', tau2) 'r^2(\tau_1) = ' sprintf('%.2f, ', rSqr1) 'r^2(\tau_2) = ' sprintf('%.2f',rSqr2)],...
                    sprintf('t_{min} = %g ms, t_{pass} = %g ms', tsMinDefo(globalCount), tsPass(globalCount))},'FontSize',8);
                ylabel('defo both')
            subplott(3,1,2);
                plot(time(maskInsideChannel),1-circFCEven(maskInsideChannel),'.-', tFine,yFitCircEven);
                ylabel('defo even');
            subplott(3,1,3);
                plot(time(maskInsideChannel),1-circFCOdd(maskInsideChannel),'.-', tFine,yFitCircOdd);
                xlabel('time [ms]');
                ylabel('defo odd');
            if flagImgExport
                print(fitFig, '-dpng',[rootdir sprintf('M%i_eval_Cell%02i_tauFit', measNum(idxMeas), sampleData(idxMeas).measData.cellData(idxCell).number)]);%, '-fillpage');
            end
        end
        if flagInterrupt
            input('next');
        end
    end
end
fprintf('\n');
%% save
if FTFITManual
    save([rootdir sprintf('M%i_eval_Cell%02i_tauFit', measNum(idxMeas), sampleData(idxMeas).measData.cellData(idxCell).number)],...
        'sampleName', 'flowRate', 'channelPosMat', 'cellNums','measNums', 'filedir', 'noOfFCs', 'rootdir', 'sampleData',...
    'rSqrsCircularityTau1','rSqrsCircularityTau2','taus1Circularity','taus2Circularity',...
    'tsMinDefo','tsPass',...
    'defoTau1Inlet','defoTau1Outlet','defoTau2Inlet','defoTau2Outlet',...
    'ak','bk','rFCOdd','rFCEven','circFCEven','circFCOdd','inertRatFCEven','inertRatFCOdd','xPos','maskInsideChannel','time','timeInsideChannel',...
    'para1', 'yFitCircEven','para2', 'yFitCircOdd', 'tFine');
else
    save([rootdir 'data.mat']', 'sampleName', 'flowRate', 'channelPosMat', 'cellNums','measNums', 'filedir', 'noOfFCs', 'rootdir', 'sampleData',...
        'rSqrsCircularityTau1','rSqrsCircularityTau2','taus1Circularity','taus2Circularity',...
        'tsMinDefo','tsPass',...
        'defoTau1Inlet','defoTau1Outlet','defoTau2Inlet','defoTau2Outlet');
end

if FTFITManual
    return;
end
%% plot results
recFig = 6;

tauRange1 = [0.1 10];
tauRange2 = [0.1 1e3];

%% tau1(Circularity) histo
recFig = recFig + 1; set(0,'CurrentFigure',recFig); clf;
maskCtrlCirc = (taus1Circularity >= 0) & (rSqrsCircularityTau1 >= 0.8);
plotHisto5( taus1Circularity(maskCtrlCirc), taus1Circularity(maskCtrlCirc), tauRange1, [sampleName ' @ ' flowRate], 'lg(tau_1) in ms', 'log10', [] );
% if flagExport
    saveas(recFig,[rootdir 'CircTau1Histo.fig']);
    print(recFig, '-dpng','-painters',[rootdir 'CircTau1Histo']);
% end
%% tau2(Circularity) histo
recFig = recFig + 1; set(0,'CurrentFigure',recFig); clf;
maskCtrlCirc = (taus2Circularity >= 0) & (rSqrsCircularityTau2 >= 0.8);%(taus2Circularity < 20) &
if length(taus2Circularity(maskCtrlCirc)) > 0
    plotHisto5( taus2Circularity(maskCtrlCirc), taus2Circularity(taus2Circularity >= 0), tauRange2, [sampleName ' @ ' flowRate], 'lg(tau_2) in ms', 'log10', [] );
end
% if flagExport
    saveas(recFig,[rootdir 'CircTau2Histo.fig']);
    print(recFig, '-dpng','-painters',[rootdir 'CircTau2Histo']);
% end
